A model colloidal fluid with competing interactions: bulk and interfacial properties 
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Using a simple mean-field density functional theory theory (DFT), we investigate the structure 
and phase behaviour of a model colloidal fluid composed of particles interacting via a pair potential 
which has a hard core of diameter a, is attractive Yukawa at intermediate separations and repulsive 
Yukawa at large separations. We analyse the form of the asymptotic decay of the bulk fluid correla- 
tion functions, comparing results from our DFT with those from the self consistent Ornstein-Zernike 
approximation (SCOZA). In both theories we find rich crossover behaviour, whereby the ultimate 
decay of correlation functions changes from monotonic to long-wavelength damped oscillatory decay 
on crossing certain lines in the phase diagram, or sometimes from oscillatory to oscillatory with 
a longer wavelength. For some choices of potential parameters we find, within the DFT, a A-line 
at which the fluid becomes unstable with respect to periodic density fluctuations. SCOZA fails to 
yield solutions for state points near such a A-line. The propensity to clustering of particles, which is 
reflected by the presence of a long wavelength 3> cr, slowly decaying oscillatory pair correlation func- 
tion, and a structure factor that exhibits a very sharp maximum at small but non zero wavenumbers, 
is enhanced in states near the A-line. We present density profiles for the planar liquid-gas interface 
and for fluids adsorbed at a planar hard wall. The presence of a nearby A-transition gives rise to 
pronounced long-wavelength oscillations in the one-body densities at both types of interface. 
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I. INTRODUCTION 

Recently colloidal systems have been synthesised in 
which the effective pair potential between the col- 
loidal particles is attractive just outside the core sep- 
aration distance, but is repulsive for larger particle 



separations 



1,2,3,4.5,6 



The long range repulsion stems from 



the colloids being charged and the short range attraction 
is generated by the depletion mechanism, 7 arising from 
the addition of non-adsorbing polymers to the solution. 8 
Such competing interactions can give rise to phase be- 
haviour that can be very different from that found in 
'simple' liquids. Theory and simulation in both two 
and three dimensions for model systems with compet- 
ing interactions predicts that such interactions can give 
rise to a state with undamped periodic density fluctu- 
ations, which indicates a transition to cluster or striped 
phases (microphase separation) ,2ii2iiiii2iiiii2ii£iIJiiiLI& In 
the cluster phase the colloids are ordered in such a way 
that there are assemblies containing tens or hundreds 
of particles and large voids between the clusters, con- 
taining hardly any particles. Similarly, in the striped 
phase in two dimensions the particles are arranged in 
parallel stripes with low density regions in-between the 



stripes 



2.14.18 



In three dimensions the stripes form a gel- 



like network of elongated clusters ^iii 7 . This behaviour 
is most striking given that the pair interactions between 
the particles are spherically symmetrical and suggests 
that such systems could be important candidates for de- 
veloping self-assembling pattern forming materials. The 
liquid-vapour phase transition may be preempted by a 
transition to a cluster or stripe phase, but if the long- 
range repulsion is not sufficiently strong, no permanent 
clusters or stripes are present, and a liquid- vapour phase 



transition is found. However, the liquid-gas coexistence 
curve is unusually flat in the critical region^ 

Whilst the bulk phase behaviour has begun to be un- 
derstood, there are no detailed studies of how the long- 
range decay of the pair correlation functions is affected 
by the competition between attraction and repulsion and 
very little is known about the properties of such fluids in 
inhomogeneous situations where the average fluid density 
is non-uniform. In the present paper we present a mean- 
field density functional theory (DFT) for a simple model 
of such systems which provides a means of analysing cor- 
relation functions and a first step in elucidating the prop- 
erties of inhomogeneous fluids composed of particles with 
competing interactions. 

The particular mode l fluid we consi der in this paper is 
that described in Refs. Il9l20l2ll22l23l in which particles 
interact via a double- Yukawa pair potential of the form: 



oo r < a 

f3v{r) — \ —eacxp(—Zi(r/a — l))/r 

+Aaexp(— Z-xirja — l))/r r > cr, 



(1) 



where (3 = 1/fcsT is the inverse temperature (we shall 
set (3 — 1), a is the diameter of the particles and the 
amplitudes e > and A > 0. For a fixed value of A, 
the parameter e plays a role somewhat akin to an inverse 
temperature. In the context of colloid-polymer mixtures 
one can envisage an athermal system at a fixed temper- 
ature and varying the chemical potential of a polymer 
reservoir thereby changing the strength of the depletion 
attraction. Indeed, in the simple Asakura-Oosawa-Vrij 
model, the depletion attraction between colloidal parti- 
cles arising from the exclusion of ideal polymers, yields 
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a well depth proportional to z p , the fugacity of the ideal 
polymers Following the spirit of this approach we 
can view our parameter e as a measure of the polymer 
fugacity. In much of the work described here we follow 
Pini et alJ-Lsl by setting the two (dimensionless) decay 
lengths Z\ = 1 and Z2 = 0.5. However, we also present 
some results for other choices of Z\ and Z2. 

A number of phase diagrams for this model system 
are displayed in Refs. and |2£j, for cases when A is 
not too large, i.e. cases when the repulsive contribu- 
tion to the pair potential v(r) is not sufficiently large 
that the liquid-gas phase separation is replaced by mi- 
crophase separation. For these cases one finds that there 
can be a substantial supercritical region in the phase di- 
agram where the compressibility is unusually large ^SP. 
In Refs. a number of theories have been 

used to investigate the bulk fluid correlation functions 
(radial distribution function g(r) and structure factor 
S(k)). Here we make a systematic study of the asymp- 
totic decay r — > 00 of g(r) in various regions of the phase 
diagram. The form of the decay of g(r) is determined 
by the poles of S(k) in the upper half of the complex k- 
planei^l We use a simple DFT, equivalent to the random 
phase approximation (RPA) of the tail potential outside 
the hardcore for the bulk fluid, to elucidate the precise 
pole structure for particles interacting with pair poten- 
tials given by Eq. Q . Due to the competing interactions 
the fluid displays a propensity to cluster formation which 
is manifest for some state points as a slowly decaying, 
long wavelength oscillatory decay of the radial distribu- 
tion function g(r) - the wavelength being related to the 
size of the clusters. Of course, at sufficiently low densi- 
ties the decay of g(r) is monotonic. We locate lines in 
the phase diagram at which the ultimate asymptotic de- 
cay of g(r) crosses over from monotonic to oscillatory. 
In addition we find, for some values of the fluid pair 
potential parameters, that there can be a crossover in 
the decay of g(r) from long wavelength oscillatory to os- 
cillatory with wavelength ~ a. We find that the pole 
structure displayed by our simple (RPA) DFT, and hence 
the crossover behaviour, is mimicked closely by the more 
accurate self consistent Ornstein Zernike approximation 
(SCOZA),^^ indicating that the simple DFT describes, 
at least qualitatively, the main features of the fluid struc- 
ture. 

For sufficiently large values of A we find that on in- 
creasing e, the peak in the structure factor calculated 
from the DFT can diverge at k — k c 7^ 0. We denote the 
line in the phase diagram at which this divergence occurs 
the 'A-line'i22i2ii2£ On the A-line the fluid becomes un- 
stable with respect to periodic density fluctuations with 
wavenumber k c , indicating that there is a phase transi- 
tion to a modulated phase (either a crystalline or striped 
phase), which may be pre-empted by a transition to a 
glassy non-ergodic stateii^iiS In mean field treatments 
the A-line encloses a region of the phase diagram where 
one would expect to find the liquid-gas critical point and 
it intersects the binodal at densities on either side of the 



critical point Mi We find such behaviour in the present 
DFT approach. 

This paper is laid out as follows: In Sec. [H] we intro- 
duce our simple DFT approach and in Sec. lIIII we describe 
briefly the implementation of the more accurate (bulk) 
SCOZA theory for determining correlation functions and 
phase behaviour. In Sec. IIVI we present results for the 
bulk phase behaviour for a range of pair potential pa- 
rameters. Sec. describes our approach for determining 
the ultimate asymptotic decay of the radial distribution 
functions and presents a number of representative results. 
In Sec. lVII we apply our DFT to calculate inhomogeneous 
fluid profiles at a planar hard wall and at the liquid-gas 
planar free interface. We conclude, in Sec. IVIII with a 
summary and discussion of our results. 



II. A MEAN FIELD DENSITY FUNCTIONAL 
THEORY 

The basis of DFT is that there exists a functional, 
f2y[p], of the fluid one body density profile p(r) such 
that the fluid equilibrium profile p(r) minimises this func- 
tional 



sn v [ P '} 



8p'(r) 



(2) 



p'(r)=p(r) 



and the minimal value of this functional is the grand po- 
tential, fl, of the fluids The grand potential functional 
fly can be written as: 



n v [p)=F[p)- / drp(r)[M-F ext (r) 



(3) 



where p is the chemical potential, V ex t(r) is a one-body 
external potential and F[p], the intrinsic Hclmholtz free 
energy functional, is a unique functional of p(r) for a 
given interaction potential. We can write T[p] — T i( i[p] + 
T ex \p\ where 



T ld [p\ = k B T / drp(r)[ln( )0 (r)A 3 ) - 1] 



(4) 



is the ideal-gas contribution, A the thermal de Broglie 
wavelength, and T ex is the (unknown) excess contribu- 
tion. Taking two functional derivatives of T ex , we obtain 
the (Ornstein-Zernike) pair direct correlation function for 
the inhomogeneous fluid:- 4 



c (2) (r,r') = -13 



5 2 T ex [p] 
Sp(r')Sp(r) 



(5) 



In the present work we approximate the excess 
Helmholtz free energy functional by 



T ex [p] = T^[p] + drV(r)p(r> 



>(|r-r'|), 



(6) 
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where f™. [p] is the reference hard-sphere Helmholtz ex- 
cess free energy functional and the contribution due to 
the remainder of the pair potential > a is treated in a 
mean-field fashion. 34 We define the 'perturbation' poten- 
tial as 



A 



(3v(r) 



r < a 
r > a. 



(7) 



Note that — e + A is the value of j3v(r) at contact, i.e. at 
r = <j + . We employ the Rosenfeld fundamental measure 
theor y 35 i 36 i 37 for F^[p]. This non-local functional, via 
(0, generates cpy (r, pb), the Percus-Yevick (PY) pair di- 
rect correlation function for a uniform hard sphere fluid. 
Thus, in bulk, Eqs. (0 and 10 together generate the fol- 
lowing simple (random phase) approximation (RPA) for 
the pair direct correlation function: 



c {2) {r;p b ) = c PY (r;pb) ~ 0v p (r), 



(8) 



where pb is the bulk fluid density. As in other RPA treat- 
ments of models with hard cores the perturbation poten- 
tial is not defined uniquely within the hard core. We 
choose the form in Eq. (JJJ) for simplicity and because 
the resulting phase diagrams are closer to those obtained 
from the SCOZA than those from other possible choices 
that we considered. With this choice the Fourier trans- 
form of c(r) = (r; pf,) can be carried out analytically 
and yields 



c(fc) - c PY (k) - (3v(k), 



(9) 



where cpy(fc), the Fourier transform of cpy(r; is 
given byt2& 



C PY (k) = - 4:71" (j 3 



a + 2/3 + 47 24 7\ • 

3 5~ Sm 1 

a + (3 + 7 2/3 + 127 247 

7t6~ 



r 



/247 _ 2/? 

with q = ka, and the coefficients 
(1 + 277) 2 



cosg 
(10) 



a 



P = 



-677(1 + ?7/2) 2 
(1 — v) 4 ' 

y(l + 2t?) 2 
2(1 -V) 4 



(11) 



depend upon rj = irpb<J 3 /6, the packing fraction, and 



0v p (k) 



4na 3 t 



Zi . 



—5- — - cosgH sing 

Zf + q 2 V 1 

4it<t 3 A ( Z-x . 



■ smq 



47T(T 3 (e - A) , 

^ (sin q — q cos q) . (12) 



Thus we have a relatively simple analytic expression for 
c(k) and therefore also for the static structure factor ob- 
tained from the Ornstein-Zernike relation^ 



S(k) = — 



1 



p b c(k) ' 



(13) 



In Fig. ^ we display S(k) calculated at a bulk fluid 
density pb<J 3 — 0.2457, the critical point density as ob- 
tained from DFT (see Sec. II V|) . for a number of dif- 
ferent values of the parameter e, for the case when 
A = 0.082, Zi = 1 and Z 2 = 0.5. We see that as e" 1 
is decreased, the structure factor develops a peak at a 
small, but non-zero, wave- vector k = k c <C 27t/ct. This 
peak indicates the propensity towards clustering in the 
fluidi 2 ! 5 ! 11 ! 14 ! 15 ! 21 ! 22 ! 23 A similar RPA approximation was 
used in Ref. ^| to account for the structure factor of a 
two dimensional fluid of particles with competing inter- 
actions close to those of the present fluid. The propensity 
towards clustering in the fluid can also be seen in Fig. [21 
where we plot the radial distribution function, g(r), for 
the same fluid at pb<r 3 = 0.2 and several values of e _1 . We 
observe the development of long wavelength oscillations 
in g(r) as e _1 is decreased. These results are obtained 
from the inverse Fourier transform of S(k). Note that 
within this route, the core condition that g(r) — for 
r < a is violated. This would not be the case if we were 
to use the test particle route to g(r) - i.e. if we treat one of 
the particles as a fixed external potential (V ex t(r) = v(r) 
in Eq. (|3J)), use the DFT to calculate the inhomogeneous 
fluid density profile, p(r), around this fixed particle, and 
divide by the bulk density to obtain g(r) = p(r)/pb- 

For sufficiently large values of A we find that on de- 
creasing e _1 at fixed pb the peak in the structure fac- 
tor that occurs at small k — k c + 1 can diverge. The 
line in the phase diagram at which this occurs is the A- 
j m6i 30 i 3i i 32 This line is shown in the inset of Fig.[S]for the 
case A = 0.082. Inside the region bounded by the A-line, 
the RPA result for S(k) is unphysical, since this becomes 
negative in a certain interval of k, and goes from negative 
to positive divergent values upon crossing the boundaries 
of the interval. 



III. THE SCOZA FOR CORRELATION 
FUNCTIONS AND THERMODYNAMICS 



The SCOZA is designed to deal with two-body poten- 
tials which, like that of Eq. Q , consist of a singular hard- 
sphere repulsive part, with diameter a, and a longer- 
ranged tail. As is customary in integral-equation theo- 
ries, this approach is based upon the Ornstein-Zernike 
(OZ) equation linking the radial distribution function 
g(r) to the pair direct correlation function c(r). A closed 
theory is obtained by supplementing the OZ equation 
with an approximate (closure) relation involving g(r) 
and c(r). In its simplest form, the SCOZA amounts to 
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setting 



28 



FIG. 1: Static structure factors S(k) for a fluid with total 
density pta 3 — 0.2457, the critical point density determined 
from the DFT (RPA), and with the parameters A = 0.082, 
Zi = 1 and Z2 = 0.5, calculated for = 10, 3, and 2 
(bottom to top). The solid lines are the results of the DFT 
(or RPA), Eqs. © - JT3J, and the dashed lines from SCOZA 
(there is no SCOZA result for e" 1 = 10). Note that as e" 1 
is decreased, a peak develops and grows at small wave vector 
k — k c 7^ 0, whereas the other peaks in S(k), determined 
primarily by correlations in the reference hard sphere fluid, 
are almost unchanged. 



1.5 - 




0.5 J 



FIG. 2: Radial distribution function, <?(r), for a fluid with 
total density p^o 3 — 0.2, with the parameters A — 0.082, 
Z\ — 1 and Z-2 = 0.5, calculated for e _1 = 5, 2, 1.7, and 1.66. 
Note that as e _1 is decreased, g(r) becomes longer ranged; 
the decay is oscillatory, with a wavelength of about 25a for 
(T x = 1.66 (a state close to the A-line for this system - see 
Fig-El ~ see the magnification in the inset. Within the present 
DFT (RPA) treatment, based on Fourier transforming 113L 
the core condition g(r) = for r < a is violated - see text. 



g(r) =0 r < a 

c(r) = K(pbi 0)v(r) r > a . 



(14) 



This closure differs from the RPA considered in Sec.HTIm 
two respects: first, the core condition on the radial distri- 
bution function is satisfied. Second, the amplitude K of 
the direct correlation function outside the repulsive core 
has not been set to K = —f3. Rather K is regarded as 
an unknown state-dependent quantity, to be determined 
in such a way that consistency between the compressibil- 
ity and the energy route to thermodynamics is enforced. 
The previous applications of SCOZA, including the study 
of fluids with competing interactions^*^ were aimed at 
thermal systems with a temperature-independent poten- 
tial, such that the phase diagram can be plotted as a 
function of density and temperature. For thermal flu- 
ids, the consistency condition amounts to requiring that 
the reduced compressibility Xrod and the excess internal 
energy per unit volume u satisfy the condition: 



d_ 

df3 



1 



Xrcd 



Pb 



dpi 



(15) 



where Xred is obtained from the compressibility sum rule, 
i.e. Xrcd = S(0), while u is obtained from the energy 
equation u = 2npl dr r 2 v(r)g(r). If the closure l|14|) 
is implemented for the correlations, the consistency con- 
dition 1(15(1 yields a closed partial differential equation 
(PDE) for the function K(pi,, 0). It should be noted that 
in the present case (3v{r) is in effect independent of tem- 
perature (see Eq. JIJ), and the phase diagrams will be 
plotted as a function of e _1 , the inverse strength of the 
attractive tail potential. For such an athermal system, 
Eq. (|15fl does not hold anymore, as a consequence of the 
trivial dependence on the temperature of the Helmholtz 
free energy, F ~ ksT. However, for any given e, the value 
of f3F of the athermal system is the same as that given by 
a temperature-independent interaction v(r) = w[3v{r) at 
a temperature fc^T = w, where f3v{r) is given in Eq. |T]l. 
and uu is an arbitrary energy scale. For this system, 
Eq. I |l 5(1 is valid. Therefore, the phase diagram and the 
correlations of the original potential QJ have been deter- 
mined by integrating the SCOZA PDE for v(r) down to 
the reduced temperature T* = ksT/w = 1 for each of 
the values of e considered. Equivalently, one may refor- 
mulate Eq. ((15(1 for the athermal system by replacing the 
inverse temperature (3 with a coupling constant which 
varies from zero to unity, whose purpose is to switch on 
the tail interaction. 

In the present case, implementing the SCOZA scheme 
is made simpler by taking advantage of the analytical 
results^ obtained for the mean spherical approximation 
(MSA) when the tail potential has a two- Yukawa form as 
in Eq. . These enable one to obtain Xred as a function 
of pi, and u. The function Xred{pb, u ) can then be used in 
Eq. 1 (15(1 by taking u instead of K as the unknown quan- 
tity. The algebraic manipulations have been described in 
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detail in Ref . |45| The resulting PDE for u supplemented 
with suitable initial and boundary conditions^ is inte- 
grated numerically. Integration of u with respect to (3 
then yields the Helmholtz free energy and hence all the 
other thermodynamic quantities. 

The SCOZA, like the RPA, yields a spinodal curve, 
i.e., a locus in the temperature-density plane where the 
compressibility diverges. As soon as one enters the region 
bounded by the spinodal, Xred is n0 longer positive, so 
that the theory ceases to be meaningful. Therefore, the 
region bounded by the spinodal must be excluded from 
the integration domain. If the temperature is below its 
critical value and one approaches the critical density p c 
either from p = or from the high-density boundary po, 
it is found that Xrcd is no longer positive for p = p\ s or 
P = P2s, where p\ s and p2 S are temperature-dependent 
densities with p% s < p c < p2 S - The integration is then 
restricted to the intervals (0, pi s — Ap) and (p2 S + Ap, po), 
where Ap is the density spacing. Within the accuracy of 
the numerical discretization, p ls and pi s give the densites 
of the spinodal curve at the temperature considered. The 
liquid-gas coexistence curve is determined by equating 
the pressures and chemical potentials on the low- and 
high-density branches of the subcritical isotherms. 



IV. RESULTS FOR PHASE DIAGRAMS 

In this section we present results for the liquid-gas co- 
existence curve (binodal) and spinodal calculated from 
both theories. Within the present DFT it is straightfor- 
ward to determine liquid-gas coexistence. By replacing 
p(r) = pi, in the Helmholtz free energy functional JfJJl, we 
obtain an expression for the bulk Helmholtz free energy 
per particle, /. The corresponding pressure can exhibit 
a van der Waals loop and we calculate the coexisting 
gas and liquid densities in the standard way. We obtain 
the spinodal from the free energy as the locus of state 
points where d 2 f/dv 2 = 0, where v = 1/pb is the specific 
volume - i.e. the locus where the isothermal compress- 
ibility diverges. This is identical to the locus of points 
where S(k — 0) diverges, where S(k) is given by the 
Ornstein-Zernike relation Eq. together with equa- 

tions lj^)l- p2l) . The self-consistency between the free en- 
ergy calculated directly from the functional and from the 
compressibility route is one of the appealing features of 
the simple theory. In Figs. we display the resulting 
DFT phase diagrams for Z\ = 1 and Z% — 0.5 and a 
number of choices of the parameter A. In Figs. 0] and [S] 
we also display the corresponding SCOZA results for the 
spinodal and binodal. We find quite good agreement be- 
tween the two theories, although the DFT being a mean 
field theory, is unable to describe correctly the shape of 
the coexistence curve in the region of the critical point. 

Within the present DFT, having calculated the spin- 
odal and binodal for a particular value of A, one can 
obtain results for all other values of A by a simple re- 
scaling of the vertical (e _1 ) axis. It is straightforward to 



show that, within the DFT, the critical density is inde- 
pendent of A. For the case A = (not displayed), the 
present theory yields a critical density within 5% of the 
value obtained from the more accurate SCOZAiiS How- 
ever, as A is increased the discrepancy increases. The 
general trend that the value of e _1 at the critical point is 
decreased as A is increased is found in both theories^ 

In Figs. and ED which refer to A = 0.082 and 0.5 
respectively, the dash dotted line denotes the A-line as 
obtained from the DFT (RPA). The A-line takes the 
shape of a loop which crosses both branches of the bin- 
odal and meets the spinodal at values of e _1 below where 
one would expect the critical point. 

When considering the SCOZA results, one should keep 
in mind that, unlike the RPA, there are regimes in which 
the theory cannot be solved. In particular, this is the 
case for values of A and e such that the competition be- 
tween attraction and repulsion is very strong, and the 
fluid is expected to form microphases. In the RPA, this 
regime is marked by the appearance of the A-line, where 
the structure factor S(k) has a singularity for k = k c 7^ 0. 
However, such a singularity is incompatible with the ful- 
filment of the core condition l|14|l if the direct correlation 
function in Fourier space c(k) is analytic on the real axis, 
which is indeed the case in SCOZA. As a consequence, 
the SCOZA fails to have solutions before a divergence in 
S(k) occurs. Because of the non-local character of the 
SCOZA PDE, this lack of convergence involves the full 
range of density. Therefore, is is not possible to locate 
the A-line using SCOZA. Another constraint is related 
to the numerical integration procedure, which requires 
the SCOZA PDE to be stable with respect to small fluc- 
tuations of the solution. However, such a requirement is 
not fulfilled when A is sufficiently large that the repulsive 
tail contribution to the potential overwhelms the attrac- 
tive one. In this case, the PDE behaves like a diffusion 
equation with a negative diffusion coefficient, and cannot 
be integrated numerically. An estimate of this intrinsic 
stability limit is provided by the condition that the spa- 
tial integral of the tail potential Air drr 2 v(r) should 
remain negative. How the above constraints affect the 
existence of the solution in the present representation, 
in which the phase diagram is studied as a function of 
the attraction amplitude e at fixed repulsion amplitude 
A, can be gleaned by considering the case A = 0.082 
reported in Fig. |3J For a temperature-independent two- 
Yukawa potential with the same inverse-range parame- 
ters Z\ = 1, Z2 = 0.5 as those considered here^ the 
liquid-gas transition is found to disappear when the rela- 
tive amplitude of the repulsion A/e is larger than 0.097. 
For A — 0.082 this requires e _1 > 1.18. When this condi- 
tion on e is met, the theory fails to converge below a cer- 
tain threshold temperature T t h before liquid-vapour sep- 
aration takes place. Whether this happens or not in the 
phase diagram of Fig. |5] depends on whether the reduced 
temperature T* — 1 lies below or above the threshold 
temperature. For A/e just above 0.097, Fig. 1 of Ref. Il9l 
shows that T£ h ~ 1.7 e which, for e _1 ~ 1.18, is indeed 



6 



4 







Fisher-Widom 






MONOTONIC 


line 








Kirkwood line (1) 






,■-*" ..'"""" " 


Kirkwood line (2) 






MONOTONIC 










: . OSCILLATORY 








^\ binodal 








spinodal " / 



0.2 0.4 3 0.6 0.8 1 

FIG. 3: Phase diagram for the case when Z\ = 1, Zi — 0.5 
and A = 0.0001 obtained from the DFT (RPA). The solid 
lines denote the binodal and spinodal. The two dotted Kirk- 
wood lines denote loci at which the asymptotic decay of h(r) 
crosses over from monotonic to oscillatory. For the present 
choices of Z\ and Z2, the Fisher-Widom line is at densi- 
ties pbu 3 > 0.9 which lie inside the liquid-solid coexistence 
region 1— 




FIG. 4: As in Fig. 01 except here A = 0.02 and we also display 
SCOZA results: the dashed lines denote the SCOZA spin- 
odal and binodal and the open circles correspond to points on 
the two SCOZA Kirkwood lines. Note that within the DFT 
(RPA) the lower Kirkwood (dotted) line is quite close to the 
spinodal which means that there is only a very small supercrit- 
ical region where the asymptotic decay of h(r) is monotonic. 
A similar scenario pertains with SCOZA - see magnification 
of the critical region in the inset. 



above unity. Therefore, the SCOZA solution disappears 
as soon as > 1.18. However, as e _1 increases, the 
reduced threshold temperature Tj* h decreases, until for 
(r 1 > 1.56 one has T£ h < 1, so that convergence for 
T* = 1 is again achieved. In the RPA picture, this cor- 
responds to being above the maximum of the A-line. For 
larger values of e _1 , the SCOZA can be solved. How- 
ever, the solution disappears again at high values of e _1 , 
when the repulsive part of the interaction dominates. For 
Z x = 1, Z 2 = 0.5, A = 0.082, the SCOZA PDE fails 
to converge for e _1 > 4.04, in fair agreement with the 
value (T l = 4.06 above which the spatial integral of the 
tail potential for r > a becomes positive. In summary, 
for A = 0.082 the SCOZA solution does not exist for 
1.18 < < 1.56 and for er 1 > 4.04. Turning now 
to the case A = 0.5 shown in Fig. [U the SCOZA solu- 
tion disappears for e _1 > 0.097/A = 0.194. However, for 
this choice of A the solution does not reappear at larger 
values of e _1 , because the region where the competition 
is strong and a A-line is present according to RPA over- 
laps with the region e _1 > 0.67 where the integral of 
the tail potential becomes positive and the SCOZA PDE 
is unstable. Therefore, for A = 0.5 only the interval 
er 1 < 0.194 is accessible to the SCOZA. We did calcu- 
late portions of the binodal and spinodal in this interval, 
finding good agreement with the DFT results. However, 
for these comparatively low values of eT 1 the relevant 
densities are either very low or very high, and are not 
displayed in Fig. 



ASYMPTOTIC DECAY OF CORRELATION 
FUNCTIONS 



In this section we determine the asymptotic decay, 
r — > 00, of the total pairwise correlation function h(r) = 
g(r) — 1 in our model. The basic procedure follows that 
in Ref. ;27. In Fourier space the Ornstein-Zernike (OZ) 
relation is given by Eq. 113|) . or equivalently 



h(k) = 



1 - pbc(k) 



(16) 



where h(k) is the three-dimensional Fourier transform 
(FT) of h(r). Inverting the FT, and noting that h(k) is 
even, we can write: 



rh(r) 



1 



1 

An 2 i 



dkke l * r h{k) 

dkke** - m (hy (17) 
, 1 - pbc(k) 



which can be evaluated by contour integration in the 
complex planeiSl We expect the singularities of h(k) to 
be simple poles. Choosing an infinite radius semi-circle 
contour in the upper half of the complex plane, we obtain 



rh(r) 



2tt ^ 



R n e 



ik n r 



(18) 
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FIG. 5: As in Figs. El and |U except here A = 0.082. The bin- 
odal and spinodal calculated from DFT (solid lines) lie close 
to the corresponding SCOZA results (dashed lines). Note 
that solutions to SCOZA do not exist for 1.18 < < 1.56 
and for e' 1 > 4.04 for this value of A - see text. Within the 
DFT (RPA) there is a A-line (dash-dotted line in magnifica- 
tion in the inset) at which the structure factor S(k) diverges 
at k = k c 7^ 0, where k c <C 2n/a. 
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FIG. 6: As in Fig. except here A = 0.5. Now the region 
of the phase diagram enclosed by the A-line (dash-dotted) 
is quite large. Furthermore, the region enclosed by the Kirk- 
wood (dotted) line is also very large - see larger scale in the in- 
set. Note that solutions to SCOZA exist only for e~ x < 0.194 
- see text. 

where R n is the residue of kh{k) for the nth pole at k = 
k n . The k n are solutions of 

1 - p b c{k n ) = 0. (19) 

In general, there are an infinite number of poles. Poles 
lying off the imaginary axis occur in conjugate pairs 
k n = ±ai + iao and such a pair contributes a damped 



oscillatory term of the form exp(— anr) cos(ai?' — 9) to 
the sum in Eq. I|18fl . Poles that lie on the imaginary 
axis, k n = ia , contribute a pure exponential term of 
the form exp(— a^r) to the sum in Eq. 118fl . The longest 
range decay of h(r) is determined by the pole(s) with the 
smallest imaginary part. If ao < an the longest range 
decay is damped oscillatory, but if an > an, then the 
asymptotic r — > oo decay of h(r) is monotonic. 

Given the analytic expression we have for c(k) within 
the DFT (RPA) treatment (Eqs. ©-©)> it is fairly 
straightforward to calculate the full set of solutions {k n } 
to Eq. (jTT?|) for the present model fluid. Thus we are 
able to determine the pole(s) with the smallest ao and 
determine the ultimate decay of h(r). 

Before presenting our numerical results, we first discuss 
a simplified model, which should provide insight into the 
origin of the leading order poles. In Eq. (7J, v p (r) takes 
a constant value for r < a. If we replace this constant 
value by the form that v(r) takes for r > a (i.e. keeping 
a double Yukawa form for all values of r) then the de- 
nominator function D(k) = 1 — pbc(k) takes the simple 
form: 

D(k) = 1 - p b c PY (k) - + ^- F , (20) 

with q = ka. We find that this simplified model dis- 
plays the same pole structure as the full version with 
v p (r) given by Eq. Q. When A = this simplifies even 
further: 

D (k) = l-p b c PY (k)-^^. (21) 

We first seek purely imaginary poles and substitute k = 
iao into Eq. Ij21|l . Provided and < 10, then cpyiicto) ~ 
cpy(0), so that: 

D (ia)*B- C (22) 

where B = [1 — pbCpy{0)] > and C = An^ept > are 
constants dependent on the state point of the fluid. There 
is only one solution to the equation Z3 (* a o) — 0, where 
D is given by When B = C/Zf then we obtain the 
solution an = 0; this is just the spinodal. Furthermore, 
the solution is bounded above, i.e. aner < Z±. 

We find that there arc an infinite number of complex 
poles. These are essentially just the poles of the hard 
sphere reference fluid. Of these, the pole with the small- 
est imaginary part has a real part ai « 2ir/a. The imag- 
inary part of this pole, ao, takes relatively large values 
at low densities and decreases as the density is increased. 
Thus we find, for the fluid with A — 0, that at low densi- 
ties the purely imaginary pole dominates the asymptotic 
decay of h(r). However, as the density is increased, one 
finds that at some point ao = a , and on further increas- 
ing the density one finds that the asymptotic decay of 
h(r) crosses over from monotonic to damped oscillatory 
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decay, with a wavelength ~ a. The line in the phase dia- 
gram defined by the locus ao — ao is termed the Fisher- 
Widom (FW) linei 2 Ti 4 ? For the case when Z\ = 1 and 
Z2 = 0.5, this line is at high densities. In Fig. [3] we dis- 
play the FW line for the case A = 0.0001 where this line 
is at densities pa 3 > 0.9 and lies within the liquid-solid 
coexistence regioniSll 

We now consider the case A / 0. For purely imaginary 
poles, k = i&o, Eq. J3U|) becomes (cf. Eq. iJHJ)): 



D{ia) « B - 



C 



E 



2^,2 



a* a. 



(23) 



where E = Ana 3 Apt, > is independent of 6tQ. When 
B = C/Zf — EjZ\ then we obtain the solution ao = 
to the equation D(iao) = 0, which corresponds to the 
spinodal. When A is sufficiently small that there is no 
A-line in the phase diagram, then we find the following 
scenario: At the critical point (or more generally, on the 
spinodal) there are two roots of D(ia,o) = 0, the so- 
lution &o — i5o = an d another at a larger value of 
ao = <5q < Zz/a. On increasing e _1 , moving away from 
the spinodal, increases in value, whilst a h decreases 
in value. As one moves further from the spinodal these 
two roots (poles) coalesce at a minimum of D(iao) in the 
interval < ao < Zija . On moving still further from the 
spinodal, we find that there are no solutions to the equa- 
tion D(iao) = 0, i.e. there are no purely imaginary poles 
in this portion of the phase diagram. Near the spinodal, 
D(icto) < for Z2/0 < ota < Z\ja. However, on mov- 
ing away from the spinodal, one finds that in this inter- 
val the function increases and its maximum touches the 
axis D(idto) = 0. On increasing e _1 further the function 
crosses the axis and there are two roots (purely imaginary 
poles) in the interval Z2/0 < ao < Z\ja. As one con- 
tinues to move away from the spinodal these two poles 
separate and very far from the spinodal one finds that 
one pole, ag — > Z^ jo and the other pole a\ — > Z^ /<J. 

Analysing the roots of D{k) — more generally, we 
find that when the two purely imaginary poles discussed 
above coalesce, they form a conjugate pair of complex 
poles at k = ±«i + iao- At the point of coalescence 
Q'i = 0, i.e. the wavelength of the oscillations is infi- 
nite. For all the relevant values of e , a\ remains small, 
such that < a\ <C 2n/a. It is this conjugate pair of 
poles which generates the oscillatory decay of h(r) with a 
wavelength 3> cr, indicating the tendency to cluster for- 
mation. The line in the phase diagram at which these 
pairs of purely imaginary poles coalesce is denoted the 
Kirkwood line^i^ 2 ^ Kirkwood^I was the first to de- 
scribe this mechanism for crossover from monotonic to 
oscillatory decay in his study of (charge) correlations in 
electrolytes. Moving away from the spinodal in the phase 
diagram, the conjugate pair of complex poles with small 
ai moves from the real axis (increasing ao). Eventually 
«i — > and the poles rejoin the imaginary axis at some 
point in the interval Z2/0 < ao < Z\/a, coalescing to 
form the second pair of purely imaginary poles described 
above. There is therefore a second Kirkwood line in the 



phase diagram. In Fig. we display the poles calculated 
numerically from Eq. H19|) with c(k) given by Eqs. ©- 
(|12[1 along a path at the critical density, pa 3 — 0.2457, 
for a fluid with A = 0.02, Z x = 1 and Z 2 = 0.5. These 
results exhibit all the features of the simpler model de- 
scribed above. 

Fig. [7Ji plots the imaginary part of the low lying poles. 
For small values of e" 1 there is a pair of purely imagi- 
nary poles (solid line in Fig. [7^) coalescing at the first 
Kirkwood point at e _1 = 2.38 and evolving as a con- 
jugate complex pair (dash dotted line) until the second 
Kirkwood point e _1 = 6.03 when a second pair of pure 
imaginary poles emerge. Fig.[7|3 shows that the real part 
of the lowest lying complex pole ot\<j remains < 0.16 
throughout the range where the asymptotic decay of h(r) 
is oscillatory. We see that the pole with oti « 2w/a, 
arising from the reference hard sphere correlations, has 



uoo 



2.8 > Z\ and therefore for all values of e 



at 



this density, this pole does not determine the asymptotic 
decay of h(r). We also display the next higher order pole 
arising from the reference hard sphere correlations, which 
has a\ « 47t/ct and ctoo « 4.3; the imaginary and real 
parts of these poles are plotted as dashed lines in Figs.[7Jt 
and 03. Note that for the poles originating from the ref- 
erence hard sphere correlations the values of ao and a\ 
change very little as e" 1 is increased. 

As the parameter A is increased, the maximal sep- 
aration between the two Kirkwood lines increases (see 
Figs. 00; one Kirkwood line moves towards the spin- 
odal (decreasing e _1 ) and the other moves away from 
the spinodal. Within the DFT the lower Kirkwood line 
meets the spinodal at a value of A ~ 0.051 (for Z\ = 1 
and Z2 = 0.5). For greater values of A a A-line appears 
in the phase diagram, enclosing the region where one ex- 
pects the critical point. The A-line corresponds to the 
situation where the imaginary part ao — > for the com- 
plex pole that has a small real component < a.\ <C 2ir/a 
at some point above the spinodal. A divergence is gen- 
erated in the static structure factor at k = k c = a±. De- 
creasing e" 1 and following the A-line, one finds that the 
value of k c decreases continuously and eventually k c — > 0, 
i.e. the A-line converts to the spinodal. For the case 
Zi = 1 and Z2 = 0.5 we could not obtain a solution 
for the SCOZA, near where one might expect to find the 
A-line, above a certain value of A. On the basis of the 
discussion at the end of Sec. IIVI the value of A above 
which SCOZA fails to converge is determined by requir- 
ing that the threshold value of the relative amplitude^ 
A/e — 0.097 is reached when the corresponding thresh- 
old temperature Tj* h ~ 1.7e is equal to unity. This gives 
A > 0.057. This bound on A is somewhat larger than 
the value A = 0.051 predicted by the present DFT. 

We confirmed that the scenarios for the pole struc- 
ture described in Fig. are also present in the SCOZA 
theory. The calculation of poles is straightforward be- 
cause in SCOZA one has analytic expressions for c(&)44 
In Figs. and we display the Kirkwood lines (open 
circles) calculated from the SCOZA alongside those from 
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FIG. 7: (a) The imaginary part, an, of the low lying poles for 
a fluid with A — 0.02, Z\ = \ and Z2 = 0.5, calculated within 
the DFT (RPA) along a line of constant density pa 3 = 0.2457, 
the critical density. The dashed lines correspond to complex 
poles arising from the reference hard-sphere correlations, the 
dot dashed line to a complex pole with a small real part 
< qi <C 2n/a and the solid lines to two purely imagi- 
nary poles. The point near e _1 = 2.38 where the first pair of 
purely imaginary poles coalesce corresponds to the first Kirk- 
wood point and the point near e _1 = 6.03 where the second 
pair coalesce corresponds to the second Kirkwood point for 
this density - see Fig. 2] (b): The real part ai for the com- 
plex poles plotted in (a). The inset is a magnification. Note 
that Q10" = at the two Kirkwood points and remains small 
in the oscillatory region. 



DFT (dotted lines) . We find good agreement between the 
results, demonstrating that our simple DFT (RPA) de- 
scribes correctly the behaviour of the poles in the model 
fluid. 

We emphasise that there are two distinct mechanisms 
for a crossover in the asymptotic decay of h(r) from 
monotonic to damped oscillatory: (i) the coalescence 
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FIG. 8: Phase diagram for the case when Z\ — 6, Z2 — 2 
and A = 0.3. The Kirkwood line (dotted) denotes the line at 
which the asymptotic decay of h(r) crosses over from mono- 
tonic to damped oscillatory with wavelength 3> a and the FW 
line the crossover from monotonic to damped oscillatory with 
wavelength ~ a. Between the two regions of oscillatory decay 
there is another (oscillatory to oscillatory) crossover line. The 
dashed line is the continuation of the Kirkwood line, although 
in this region of the phase diagram it denotes crossover be- 
tween higher order poles, i.e. those with larger values of an 
than the leading order pole. There is also a second Kirkwood 
line located just above the critical point so that in the imme- 
diate vicinity of the critical point h(r) decays monotonically 
(OZ-like). However, on the scale of this figure, this second 
Kirkwood line is indistinguishable from the spinodal. 



of two purely imaginary poles to form a conjugate pair 
of complex poles defines a point on the Kirkwood line 
and (ii) the simpler Fisher-Widom (FW) mechanism, 
whereby on decreasing the fluid density a purely imagi- 
nary pole descends down the imaginary axis and acquires 
an imaginary part do smaller than that of the dominant 
complex pole with imaginary part ao, so that the ulti- 
mate decay becomes monotonic. As mentioned earlier, 
when Z\ = 1 and Z 2 = 0.5, the FW line is at high densi- 
ties, inside the region where we expect the solid phase to 
be the equilibrium state - see Fig.0 However, for larger 
values of Z\ and we find that the FW line moves 
to lower densities, where it can intersect the Kirkwood 
lines. In Fig. [S] we display the phase diagram for a fluid 
with Z x = 6, Z 2 = 2 and A = 0.3. We find that when 
the Kirkwood and FW lines intersect, another kind of 
structural crossover line appears in the phase diagram: 
there is a crossover from damped oscillatory decay with 
one wavelength to damped oscillatory decay with another 
wavelength4§*21 In the present case the crossover is from 
decay with a long wavelength of many particle diame- 
ters (clustering in the fluid) to a decay with a wave- 
length ~ a as e _1 is increased, see Figs. and El Note 
that for the choice of parameters corresponding to the 
phase diagram in Fig. EI the liquid-gas transition may be 
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weakly metastable. However, this does not affect the de- 
cay crossover behaviour, which is in a region of the phase 
diagram far from any phase transitions^ 

The pole analysis described above points to those re- 
gions of the phase diagram where clustering might oc- 
cur, i.e. the region between the Kirkwood lines. How- 
ever, so far we have determined only the form of the 
asymptotic decay of h(r): whether it is monotonic or os- 
cillatory as r — > oo. One must also calculate the ampli- 
tude of the longest wavelength oscillatory contribution to 
h(r) in order to assess how significant clustering is in the 
fluid. This is straightforward within the present DFT. 
The amplitude is determined by the residue of the pole 
(see Eq. (JT^J), which is given by: 



k n c(k n ) 
Pbc'{k n ) 



(24) 



where the prime denotes the derivative with respect to 
k. For states just below the upper Kirkwood line the 
amplitude is generally quite small. For example, when 
A = 0.02, for the fluid with Z x = 1 and Z 2 = 0.5 
(see phase diagram in Fig. @}, at the state point with 
Pbcr 3 — 0.25 and (T 1 = 5, there is a complex pole 
with ol\<7 = 0.0909 and olqu — 0.613. This makes 
a damped oscillatory contribution to h(r) of the form 
^4exp(— a r) cos(air — 9)/r with amplitude A — 0.0994 
and phase = 0.992, i.e. this term gives quite a small con- 
tribution to h(r), despite being the term that determines 
the ultimate asymptotic decay. However, further below 
the Kirkwood line the amplitude A can become larger 
and in the vicinity of a A-line, A can be very large so that 
the pole provides a large contribution to g(r) = 1 + h(r) 
- see for example the results in Fig. [21 Indeed, when the 
A-line is reached, the amplitude A diverges since the pole 
k n in Eq. Ij24(l corresponds to a maximum of c(k) on the 
real axis and c'(fc„) vanishes. We will return to this point 
in Sec.lvTTl 



FIG. 9: The imaginary part, ao, of the lowest lying poles for a 
fluid with Zi = 6, Z2 — 2 and A — 0.3, calculated within the 
DFT (RPA) along lines of constant density: (a) pa 3 — 0.15, 

(b) the critical density pa 3 = 0.2457 and (c) pa 3 = 0.4. The 
dashed lines correspond to complex poles arising from the 
reference hard-sphere correlations, the dot dashed line to the 
complex pole with a small real part < ai -C 27r/a, and 
the two solid lines to two purely imaginary poles. The point 
where the two lowest lying pure imaginary poles coalesce cor- 
responds to a Kirkwood point - see (a). In (b) this coales- 
cence occurs between higher order poles. The crossover at 
(T x = 1.79 is between complex poles with different oscillatory 
wavelengths - see the line in Fig. |H| The subsequent crossover 
at e _1 = 2.42 is from oscillatory to monotonic decay. In 

(c) there is oscillatory to oscillatory crossover at (T x = 0.75. 
The pole with a\o ~ 2n remains the dominant one for larger 
e _1 . Note that in (a)-(c) there is a second Kirkwood point at 
(T l ~ 0.4, very near to the critical value of e _1 , so that at 
smaller values of e~ (close to the spinodal) there are two 
purely imaginary poles, similar to the results displayed in 
Fig. EJi- 



VI. INHOMOGENEOUS FLUID: DENSITY 
PROFILES AT INTERFACES 

In this section we turn attention away from the uni- 
form fluid to situations where the one-body density is 
spatially varying. We consider two separate cases: a) the 
fluid adsorbed at a hard wall and b) the liquid-gas inter- 
face. In both cases the pole analysis for the decay of h(r) 
is directly relevant to the decay of inhomogeneous fluid 
one-body density profiles; the same pole(s) which deter- 
mine the asymptotic decay of h(r), determine the asymp- 
totic decay of the one-body density profiles i2i For den- 
sity profiles which vary only in one (Cartesian) direction, 
e.g. when the external potential V ex t(r) = V ex t(z), then 
provided V ex t(z) is sufficiently short ranged the longest 
range decay of the profile into bulk at z = 00 takes the 
form: 

p(z) - p b — p b A w exp(-a z), z — > 00 (25) 
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FIG. 10: Density profiles p(z) at a planar hard wall for a fluid 
with parameters A = 0.082, Z\ = 1 and Zi — 0.5 - see phase 
diagram in Fig. |S] The profiles are calculated for a fixed bulk 
density pba 3 = 0.2, for = 1.7, 2, 3 and 16. In the inset 
we plot In \p{z)rj i — pb<r 3 \. For = 16 the asymptotic decay 
of p(z) is monotonic (exponential). The other three states lie 
on the oscillatory side of the Kirkwood line. Note that the 
inverse decay length increases as e _1 is reduced and the A- line 
is approached. For a given value of e _1 the wavelength of the 
oscillations is the same as that of g(r) in Fig. 
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FIG. 11: Density profiles for the free gas-liquid interface cal- 
culated for various values of e _1 , for a fluid with parameters 
A — 0.082, Zi — 1 and Z2 = 0.5 - see phase diagram in 
Fig. |S] This fluid has a A-line close to where one expects 
the critical point. We find that as e _1 is increased the den- 
sity difference between the coexisting phases decreases and 
the profiles become more structured, i.e. the oscillations be- 
come more pronounced on both sides of the interface. This is 
opposite to what occurs in simple fluids and reflects the the 
proximity of the coexisting states to the A-line. 



when the pole with the smallest value of ao is purely 
imaginary and 

p(z) — pb ~ pbA w exp(— aoz) cos(a\z — 6 W ), z — > 00 

(26) 

when the pole with the smallest ao is complex. The 
poles k n = ictQ or k n = ±ai +icto correspond to the bulk 
density pb^ 

Using the DFT described in Sec. [HJ we have deter- 
mined the one-body density profiles p(z) for the two 
cases mentioned above. Minimizing our approximation 
for fly [p] yields an Euler-Lagrange equation for the equi- 
librium density profile p(z) which can be solved by means 
of a standard iterative scheme. In Fig.^H we display the 
one-body density profiles at a planar hard wall with 



V ex t{z) 



OO Z < a/2 

z > a/2. 



(27) 



The fluid has parameters A = 0.082, Z 1 = 1 and Z 2 = 0.5 
(see phase diagram in Fig. 15)1. The density profiles are 
calculated for a fixed bulk density p^cr 3 = 0.2, for e _1 = 
16, 3, 2 and 1.7. For e^ 1 = 16 the state point is above the 
Kirkwood line, and the bulk asymptotic decay of h(r) is 
monotonic of the form in Eq. (1251) . This is confirmed in 



the inset of Fig. 1101 which plots In \p(z)a 3 — /9&o" 3 | versus 
z; when the decay of the density profile is of the form 
in Eq. 112511 the decay of the profile in this representation 
is a straight line. The other profiles are calculated for 
state points inside the Kirkwood line, with increasing 
proximity to the A-line. From the inset to Fig. 1101 we see 
that the decay of these profiles is indeed of the form in 
Eq. (|26|l . The profile for e _1 = 1.7, a state point near to 
the A-line, has a small decay length and oscillates with 
a long wavelength of about 23a. The density profiles for 
state points in the vicinity of the A-line are significantly 
different from those far away. These profiles should be 
compared with the radial distribution functions g(r) in 
Fig. [3 where the same trends with e _1 are observed. 

It is striking that as e -1 is reduced the contact density 
p(a/2) reduces dramatically and for = 1.7 there is a 
region where the density near the wall is significantly de- 
pleted below bulk out to distances of about 4a. That the 
density at contact should become smaller follows from 
the hard- wall sum rule: ksT ' p(a /2) = p, where p is the 
pressure of the bulk fluid far from the wall. Decreas- 
ing (T 1 corresponds to turning on more of the attractive 
interaction thereby decreasing p. 

In Fig. ^] we display the density profiles for the free 
planar gas-liquid interface of the same fluid. Profiles 
are calculated for coexisting states with e^ 1 = 1.4, 1.5, 
1.6 and 1.64. This fluid exhibits a A-line enclosing the 
region of the phase diagram where one expects the crit- 
ical point, see the inset to Fig. [3] States outside the 
A-line, i.e. for e _1 < 1.65, correspond to conventional 
coexistence between (disordered) liquid and gas. How- 
ever, as the four coexisting states that we consider lie 
inside the Kirkwood line, within our mean field DFT 
we expect the asymptotic decay of the density profiles 
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to be damped oscillatory on both sides of the interface. 
This is indeed what we find. Such behaviour was found 
previously in DFT studies of fluid-fluid interfaces of the 
binary Gaussian core model4& What is striking about 
the present results and what makes them different from 
results for other models is that as e _1 is increased and 
the density difference between the coexisting phases de- 
creases, the profiles become more structured at the in- 
terface. Normally, the interfacial density profile becomes 
less structured when the difference in density between 
the coexisting phases decreases. The difference is due to 
the fact that the present system exhibits a A-line and the 
coexisting densities become closer to this line as e _1 is 
increased - see Fig. |5j The proximity of the A-line im- 
plies long-wavelength oscillatory density profiles decaying 
slowly into bulk on either side of the interface. It should 
also imply that the propensity to cluster formation in 
both phases is stronger. Clustering should manifest it- 
self in highly structured interfacial profiles. 



VII. DISCUSSION 

Using a simple DFT and the SCOZA integral equation 
theory we have investigated the bulk structure and phase 
behaviour of a model colloidal fluid with competing in- 
teractions. In particular, we have examined in detail the 
asymptotic decay of the radial distribution function g(r) 
and find a rich variety of decay scenarios. The presence 
of clustering in the fluid, which occurs for sufficiently 
large values of the repulsion amplitude A in the pair po- 
tential Eq. CJ, i s reflected in a long wavelength (3> a) 
damped oscillatory decay of g(r). For small values of 
A the region of the phase diagram in which the decay 
of g(r) is long wavelength oscillatory is bounded above 
and below by two Kirkwood lines - see Figs.Oland^ On 
these lines the asymptotic decay of g(r) crosses over from 
damped oscillatory to monotonic (exponential) decay. As 
A is increased we find within the mean-field DFT (RPA) 
that the lower Kirkwood line is replaced by a A-line (see 
Figs. [S] and EJ, indicating a transition to a phase with un- 
damped periodic density fluctuations - e.g. a cluster or a 
stripe phase. Note that the appearance of a region in the 
phase diagram with long wavelength oscillatory decay of 
correlation functions is not restricted to fluids with com- 
peting interactions. Such behaviour is to be expected in 
any system with competing interactions!^^ In Ref. 
the authors investigate the form of the correlation func- 
tions in frustrated 0(n) spin systems. They find that as 
the strength of the frustration parameter (the analogue 
of the parameter A in the present system) is increased, 
a region in the phase diagram opens up in which the de- 
cay of spin-spin correlations is oscillatory, crossing over 
above and below to regions of monotonic decay, in much 
the same way as in the present fluid system. 

In general, the results of our simple mean field DFT 
(RPA) are in good qualitative and sometimes quantita- 
tive agreement with the more sophisticated (and more 



accurate) SCOZA theory. This reflects the fact that 
our focus has been on state points where the parame- 
ters A, e < 1. Recall that (— e + A)ksT is the value of 
the pair potential at contact, r = cr + , so that we have 
focused on cases where the portion of the pair poten- 
tial that we treat in mean field fashion has an amplitude 
< ksT . We expect the present mean field DFT to be 
less reliable when (e — A) > 1, or when the width of the 
attractive portion of the pair potential becomes narrow, 
i.e. when Z\ becomes very large. In the present study we 
have, in the main, avoided these regimes; an exception is 
the system described in Fig. [8] where Z\ = 6. 

In Sec. lVII we presented results from DFT for inhomo- 
geneous one-body fluid density profiles. The oscillatory 
density profiles in Figs. and ^2 correspond to state 
points between the Kirkwood and the A-line in the phase 
diagram of Fig. [3] where the system is a disordered fluid, 
albeit with a propensity towards cluster formation. For 
state points in the vicinity of the A-line we find some 
very striking density profiles with pronounced long wave- 
length (» a) oscillations. It is important to enquire 
whether such behaviour pertains beyond the mean-field 
(MF) treatment we present here. 

Determining the existence or non-existence of a A-line 
is a problem that arises for other types of fluids, e.g. mi- 
croemulsions near the transition to a lamellar phase and 
smectic liquid crystals. Recently much effort has focused 
on a possible A-transition for charge ordering in the re- 
stricted primitive model (RPM) in which the two species 
of ions are modelled by equal sized hard spheres carrying 
equal and opposite charges. 30 It is well known that in 
the RPM the location of the A-line in the phase diagram 
depends sensitively on the choice of (MF) approxima- 
tion - see Ref. |30| and references therein. For example, a 
simple DFT yields a A-line of continuous transitions be- 
tween a uniform disordered phase and a charge ordered 
phase, whereas the more sophisticated MSA predicts no 
A-transition. For the former case Ciach et alm^ studied 
correlation functions in some detail and found that on 
approaching the A-line, from the disordered side, both 
the decay length and the amplitude of the charge-charge 
pair correlation function diverge as t -1 / 2 , where r mea- 
sures the 'distance' of the state point from the A-linc. 
The derivation of Ciach et al. is based on a particular 
form for the (charge-charge) inverse correlation function. 
However, their results generalise straight forwardlji^S to 
any fluid system where the dominant poles approach the 
real axis at Re(k n ) = k c ^ and the other poles are well 
removed. As this situation pertains in the present DFT 
(RPA) treatment of our model fluid we expect that close 
to the A-transition 

rh[r) ~ ^4exp(— a^r) sin(fc c r), r — > oo (28) 

with ao ~ t 1 / 2 and A ~ t -1 / 2 , where again r measures 
the distance from the A-line. This result arises also in a 
simplified treatment based on the small- fc expansion of 
the direct correlation function: l/S(k) = 1 — pc(k) ~ 
a + bk 2 + ck 4 . For a > 0, c > 0, the condition that 
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l/S(k) be positive definite for every fc > is violated 
for b < 0, b 2 — Aac > and the appearance of the A-line 
corresponds to the marginal case b 2 — Aac = 0. As the 
A-line is approached from the disordered side, one has 
b < 0, b 2 — Aac — > 0~, which via contour integration, 
gives Eq. (2BJ with A ~ 1/y/Aac- b 2 . Clearly the MF 
treatment implies the unphysical result that pair corre- 
lations are unbounded on approaching the A-line. What 
occurs beyond MF theory? For the RPM a field theo- 
retic treatment^ based on the approach of Brazovskii, 52 
shows that incorporating fluctuations leads to the disap- 
pearance of the A-transition. Rather a first order tran- 
sition to a charge-ordered (crystalline) phase occurs at a 
temperature below that of the original MF A-transition. 
(Note that the A-transition is absent in simulation stud- 
ies of the continuum RPM.) Moreover the charge-charge 
correlation function changes smoothly near the original 
MF A-transition; the decay length and amplitude vary 
continuously^ 

One might infer that the very striking long-ranged de- 
cay of g(r) observed in Fig. [21 for eT l = 1.66, at a state 
point rather close to the A-line in Fig. [S] is an arti- 
fact of the MF treatment since including fluctuation ef- 
fects would remove the A-transition, replacing this with 
a first-order transition to some ordered state (cluster or 
stripe phase, perhaps). For state points somewhat fur- 
ther from the A-line, such as = 2.0, the long wave- 
length oscillations still develop in g{r) but these are more 
highly damped. Nevertheless, the corresponding struc- 
ture factor S(k) exhibits a pronounced peak at a non- 
zero wavenumber k c - see Fig. Q] The fact that SCOZA 
yields a very similar structure factor gives us some con- 
fidence in the results of the mean field DFT (RPA) for 
such state points. On the other hand, unlike the RPA, 
SCOZA does not yield a A-line, rather it fails to converge 
in the regime in which a A-line is expected according to 
the RPA. As the convergence limit is approached, the 
cluster peak in S(k) may reach values much larger than 
unity, but nevertheless the peak height does not grow ar- 
bitrarily, since a singularity in S(k) at nonvanishing k is 
incompatible with the core condition Eq. I|14f> . which is 
fulfilled in SCOZA. This can be regarded as an indication 
that the A-transition will be removed when one goes be- 
yond the MF approximation. It is tempting to introduce 
a rough and ready criterion, analogous to the Hansen- 
Verlet criterion for freezing, that establishes when we ex- 
pect the MF treatment to fail or to imply the onset of 
ordering. Such a criterion might be when the peak height 
S(k c ) ~ O(l0 r ). 

Returning to the one-body profiles, we note that none 
of the bulk statepoints in Fig. ^] for adsorption at the 



planar hard wall are particularly close to the A-line and 
we might expect our DFT results to be at least qualita- 
tively correct. However, in Fig.^Jfor the density profiles 
at the planar liquid-gas interface it is clear that the strik- 
ing increase in structuring as is increased is a direct 
consequence of the close proximity of the A-line. It is 
unlikely that such behaviour could survive beyond MF. 
A possible scenario when fluctuations are included is the 
phase diagram of Fig. 0] where there are two Kirkwood 
lines and the critical point lies below the lower line (2). 
The MF DFT would yield monotonically decaying pro- 
files on both sides of the interface for coexisting states 
just below the critical point but for 1.9 < < 2.25 one 
would expect damped oscillatory profiles on both the liq- 
uid and the gas side. For states with e _1 < 1.9, where the 
upper Kirkwood line (1) intersects the gas binodal, the 
asymptotic decay on the gas side would be monotonic. 
Of course, the MF DFT omits effects of capillary wave 
like fluctuations that act to erode the amplitude of the 
oscillations; this mechanism is discussed in Ref.|2(|for the 
case of model colloid-polymer mixtures. We believe that 
much more work is required to understand the nature of 
the fluid correlations in such systems, particularly in re- 
gions of the phase diagram where the present MF theory 
predicts a A-line. Indeed, some of us are already en- 
gaged in such studies using Monte-Carlo simulations and 
various integral equation theories. Of course, one could 
also employ DFT to investigate possible spontaneous or- 
dering, i.e. whether there are solutions corresponding to 
non-uniform (e.g. periodic; stripe or cluster) phases that 
have a lower free energy than those corresponding to a 
uniform phase. 

We conclude by speculating that since colloidal sys- 
tems with competing interactions display modulated 
structures, these systems will have interesting optical 
properties and may have applications in display technolo- 
gies. 
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